About

For this project we will be looking at LiDAR-generated Canopy Height Models and looking at uncertainty compared to field measurements of plot-level and tree-level characteristics.

For this project, I am going to be working with data from NEON Domain 01, Harvard Forest Field Site (HARV). I downloaded my data from NEON’s figshare site.

Specific characteristics we will be looking at include… *Fill this part in

First, I’m going to load the libraries I need to work with the data.

Also, added additional setup steps in this code chunk.

library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# also set stringsasFactors to false
options(stringsAsFactors = FALSE)

Next, I’m going to load in the canopy height model

site_chm <- raster("../NEONdata/D01-Massachusetts/HARV/2014/lidar/HARV_lidarCHM.tif")

# take a look at it
plot(site_chm, main="Canopy Height Model\nHARV Field Site")

hist(site_chm, main="CHM Histogram\nHARV Field Site")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 10% of the raster cells were used. 100000 values used.

# set 0 values to NA
site_chm[site_chm == 0] <- NA

# take another look at the histogram
hist(site_chm, main="CHM Histogram\nHARV Field Site",
     xlab="Height (m)")

Then, I’m going to load in the insitu data

site_insitu <- read.csv("../NEONdata/D01-Massachusetts/HARV/2014/insitu/veg_structure/D01_2012_HARV_vegStr.csv", 
                        header=TRUE, sep=",")

# let's make sure the data came in okay
head(site_insitu)  # lots of NA values
##   site_id       sitename plotid easting northing taxonid
## 1    HARV Harvard Forest    A01      NA       NA      NA
## 2    HARV Harvard Forest    A01      NA       NA      NA
## 3    HARV Harvard Forest    A01      NA       NA      NA
## 4    HARV Harvard Forest    A01      NA       NA      NA
## 5    HARV Harvard Forest    A01      NA       NA      NA
## 6    HARV Harvard Forest    A01      NA       NA      NA
##         scientificname individual_id pointid individualdistance
## 1 Vaccinium corymbosum            NA      41                 NA
## 2 Vaccinium corymbosum            NA      42                 NA
## 3          Acer rubrum            NA      43                 NA
## 4     Castanea dentata            NA      44                 NA
## 5     Castanea dentata            NA      45                 NA
## 6     Tsuga canadensis            NA      46                 NA
##   individualazimuth dbh dbhheight basalcanopydiam basalcanopydiam_90deg
## 1                NA 0.5        NA              NA                    NA
## 2                NA 0.4        NA              NA                    NA
## 3                NA 5.3        NA              NA                    NA
## 4                NA 6.5        NA              NA                    NA
## 5                NA 8.3        NA              NA                    NA
## 6                NA 7.5        NA              NA                    NA
##   maxcanopydiam canopydiam_90deg canopy_diameter vertical_position
## 1            NA              0.0             0.7        Understory
## 2            NA              0.0             0.8        Understory
## 3            NA              5.4             3.7        Understory
## 4            NA              0.0             5.0        Understory
## 5            NA              7.7             0.0        Understory
## 6            NA              0.0             0.0        Understory
##   stem_height stemremarks stemstatus canopyform livingcanopy inplotcanopy
## 1         1.5          NA         NA         NA           NA           NA
## 2         1.6          NA         NA         NA           NA           NA
## 3        10.5          NA         NA         NA           NA           NA
## 4         7.0          NA         NA         NA           NA           NA
## 5         2.4          NA         NA         NA           NA           NA
## 6         4.4          NA         NA         NA           NA           NA
##   materialsampleid dbhqf stemmapqf plant_group common_name aop_plot
## 1               NA    NA        NA          NA          NA       NA
## 2               NA    NA        NA          NA          NA       NA
## 3               NA    NA        NA          NA          NA       NA
## 4               NA    NA        NA          NA          NA       NA
## 5               NA    NA        NA          NA          NA       NA
## 6               NA    NA        NA          NA          NA       NA
##   unique_id
## 1        NA
## 2        NA
## 3        NA
## 4        NA
## 5        NA
## 6        NA
# look at a histogram of the heights
hist(site_insitu$stem_height, main="In Situ Histogram\nHARV Field Site",
     xlab="Height (m)")

In the insitu data, it looks like there’s a lot more shorter vegetation than was captured in the CHM.

Load in the plot locations

# loading using readOGR

site_plots <- readOGR("../NEONdata/D01-Massachusetts/HARV/vector_data/",
                      "HARV_PlotCentroids")
## OGR data source with driver: ESRI Shapefile 
## Source: "../NEONdata/D01-Massachusetts/HARV/vector_data/", layer: "HARV_PlotCentroids"
## with 21 features
## It has 25 fields
# look at the plot locations
plot(site_plots)

# overlay plot locations on CHM
plot(site_chm, main="CHM with In Situ Locations\nHARV Field Site")
plot(site_plots, add=TRUE)

Following along with NEON tutorial on creating field data polygons from centroids.

Add link here

Set up plot dimensions

I’m going to set up square plots around the plot centroid. This process involves setting a radius, which will be half the length of one side of the square.

This is based on the assumption that plots are oriented north and not rotated.

# first set radius

radius <- 20  # in meters

# then define the plot boundaries using the radius
yPlus <- site_plots$Y + radius
xPlus <- site_plots$X + radius
yMinus <- site_plots$Y - radius
xMinus <- site_plots$X - radius

Extract Plot IDs from centroids

ID <- site_plots$plotID

# let's look at our plot IDs
ID
##  [1] "HARV_015" "HARV_033" "HARV_034" "HARV_035" "HARV_036" "HARV_037"
##  [7] "HARV_038" "HARV_039" "HARV_040" "HARV_041" "HARV_042" "HARV_043"
## [13] "HARV_044" "HARV_045" "HARV_046" "HARV_047" "HARV_048" "HARV_049"
## [19] "HARV_050" "HARV_051" "HARV_052"

Now I’m going to calculate square polygon coordinates for each plot centroid

# here I am just making a list of the coordinates of each corner of the square
# will follow x, y, x, y, etc.
# later we will turn this into the 4 corners of the square the bounds the plot
# this requires 5 corners to draw the square - why are there 6?

square <- cbind(xMinus, yPlus, 
                xPlus, yPlus, 
                xPlus, yMinus, 
                xMinus, yMinus, 
                xMinus, yPlus,
                xMinus, yPlus)

Create spatial polygons using the coordinates

# Create a function to do this
polys <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)  # take a list and create a matrix
  Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID),proj4string=CRS(as.character("+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")))

Create shapefile

This is an important step to extract the data.

polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

Plot this with our CHM

plot(site_chm, main="CHM with square in situ plots\nHARV Field Site")
plot(polys.df, add=TRUE)

Use square buffers to extract information from LiDAR data

# extracting max height from CHM at field data sites
chm_maxheight <- extract(site_chm, polys.df,
                       fun=max, sp=TRUE,  # sp for return spatial object
                       stringsasFactors=FALSE)

# extract mean height from CHM at field data sites
chm_meanheight <- extract(site_chm, polys.df,
                       fun=mean, sp=TRUE,  # sp for return spatial object
                       stringsasFactors=FALSE)

# check outputs
chm_maxheight@data
##                id HARV_lidarCHM
## HARV_015 HARV_015         19.73
## HARV_033 HARV_033         30.06
## HARV_034 HARV_034         29.86
## HARV_035 HARV_035         23.20
## HARV_036 HARV_036         29.09
## HARV_037 HARV_037         25.58
## HARV_038 HARV_038         28.87
## HARV_039 HARV_039         25.74
## HARV_040 HARV_040            NA
## HARV_041 HARV_041         27.34
## HARV_042 HARV_042         26.47
## HARV_043 HARV_043         26.24
## HARV_044 HARV_044         23.92
## HARV_045 HARV_045         25.82
## HARV_046 HARV_046         20.02
## HARV_047 HARV_047         22.85
## HARV_048 HARV_048         26.85
## HARV_049 HARV_049         23.51
## HARV_050 HARV_050         28.59
## HARV_051 HARV_051         24.97
## HARV_052 HARV_052         28.87
chm_meanheight@data
##                id HARV_lidarCHM
## HARV_015 HARV_015      15.76019
## HARV_033 HARV_033      24.50129
## HARV_034 HARV_034      22.06391
## HARV_035 HARV_035      16.70113
## HARV_036 HARV_036      24.13921
## HARV_037 HARV_037      19.59647
## HARV_038 HARV_038      21.68227
## HARV_039 HARV_039      21.75469
## HARV_040 HARV_040            NA
## HARV_041 HARV_041      19.79854
## HARV_042 HARV_042      21.94006
## HARV_043 HARV_043      19.49369
## HARV_044 HARV_044      17.43806
## HARV_045 HARV_045      21.03539
## HARV_046 HARV_046      12.43345
## HARV_047 HARV_047      15.44984
## HARV_048 HARV_048      19.45156
## HARV_049 HARV_049      18.45651
## HARV_050 HARV_050      20.35642
## HARV_051 HARV_051      20.80787
## HARV_052 HARV_052      19.81688
unique(site_insitu$plotid)
##  [1] "A01" "A04" "A05" "A06" "A07" "A08" "A09" "A12" "A13" "A16" "A17"
## [12] "A18" "A19" "A20" "A31" "A34" "A10" "A15"

Extract information from insitu data

# extract max ht by plot from insitu data with dplyr
insitu_maxheight <- site_insitu %>% 
  group_by(plotid) %>%
  summarise(max(stem_height))

names(insitu_maxheight) <- c("plotid","max_height")

hist(insitu_maxheight$max_height)
hist(chm_maxheight@data$HARV_lidarCHM)

# quick look at data summary
summary(fieldmaxht.sjer)
summary(height.sjer@data$SJER_lidarCHM)

# update names
names(fieldmaxht.sjer) <- c("plotid","field_maxht")